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ABSTRACT 


This paper discusses further improvements made to advance the current Boundary Layer Integral Matrix 
Procedure - Version J (BLIMPJ) containing previously modeled simplified calculation methods by accounting 
for condensed phase, thick boundary layer and free-stream turbulence effects. The condensed phase effects 
were included through species composition effect considered via input to the code and through particle 
damping effect considered via a turbulence model. The thrust loss calculation procedure for thick boundary 
layer effects was improved and the optimization of net thrust with respect to nozzle length was performed. 
The effects of free-stream turbulence were approximately modeled in the turbulence model. 


NOMENCLATURE 


e 

= 

Kinematic Eddy Viscosity, ft 2 /sec 

v f = 

Volume of each particle, ft 3 

P 

= 

Coefficient of Viscosity, lbm/(ft-sec) 



MF 

= 

Mass Fraction 

Subscript 


n 

= 

Number Density, no./ft 3 

D = 

Diffusivity 

V 

= 

Kinematic Viscosity, ft 2 /sec 

e = 

Edge 

P 

= 

Density, lbm/ft 3 

/ = 

Fluid 

T 


Shear Stress, lbf/ft 2 

mix = 

Mixture 

It 

= 

X Component of Velocity, ft /sec 

P = 

Particle 

V 

= 

Y Component of Velocity, ft /sec 



W 

= 

Mass Flux, lbm/(ft 2 -sec) 




INTRODUCTION 


BLIMPJ [1] has been identified by the propulsion community as the rigorous boundary layer program in 
connection with the existing J ANNAF reference programs such as ODE, ODK and TDK-BLM (all in Ref. [2]). 
The improvements made in Refs. [3], [4] and [5] have potential applications in the design of the future Orbit 
Transfer Vehicle (OTV) engines. These engines will utilize a high chamber pressure expander cycle operation 
mode which primarily depends on the heat energy transmitted from the combustion products through the 
thrust chamber wall. The larger the regenerative heat transfer, the higher the chamber pressure, which, in 
turn, permits larger area ratio nozzles. The heat transfer to the nozzle wall is affected by such variables as wall 
roughness, relaminarization and the presence of particles in the boundary layer flow. The motor performance 
loss calculation for these nozzles with thick boundary layers is inaccurate using the conventional JANNAF 
procedure. Thus, engineering procedures are needed to model these effects adequately. 
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Some of the results covering the above areas were presented at the JANNAF Combustion Meeting held 
in Monterey in 1987 [6], The most recent effort concentrated on modifying the remainder of the previously 
modeled methods. Some of these required modification of the turbulence models currently coded in BLIMP J. 
Many details have been removed in this paper for conciseness , but the reader is urged to refer to Ref. [5] for 
details. This paper is divided into three sections which describe simplified analytical formulations to include 
the following effects: 

• Presence of condensed phase (or particles) in the boundary layer 

• Thrust decrement calculation procedure for a thick boundary layer 

• Free-stream turbulence in a thrust chamber 

These effects are described in the following sections along with a set of recommendations. 


EFFECTS OF CONDENSED PHASE 


BACKGROUND 

Since the OTV nozzles will be designed to operate at very low ambient pressures and flow will expand 
very rapidly in the nozzle, there is every likelihood of forming a two-phase fluid containing water droplets 
and ice particles in the nozzle. The effect of particle loading on skin-friction and heat transfer characteristics 
at the wall is the subject investigated in this task. It has been found from studies on these high area ratio 
nozzles, which are regeneratively cooled for part of the nozzle and radiation-cooled for the last portion of the 
nozzle, that heat transfer decreases drastically downstream of the throat and is almost negligible a moderate 
distance aft of the throat. So, the formation of water droplets and ice particles occurring in this part of the 
nozzle might not greatly impact the magnitude of heat transfer and skin friction. However, the presence 
of particles in the boundary layer would change the displacement and momentum thicknesses and, in turn, 
would alter the boundary layer losses and ultimately affect performance of the thrust chamber. An analytical 
approach containing a modification of the turbulence model is given below to quantify the particle effects in 
the boundary layer. 

TURBULENCE MODEL REPRESENTATION 

Heat transfer will be affected at the wall due to the presence of particles in the boundary layer because 
of the following two significant effects. 

1. The eddy viscosity in the boundary layer is expected to decrease below its clean flow value because 
of particle damping in the boundary layer. 

2. The specific heat of the two-phase flow will be modified because of the loading, resulting in a change 
in heat transfer to the wall. 

For small particles, it has been shown by Soo and Tien [7] that the particle diffusivity is of a similar 
order as the eddy diffusivity of the stream. Thus, the turbulence model needs no modification. However, the 
remaining particle loading effect which was quantified and considered by Tien [8], and Farbar and Morley [9] 
has been reported by Praharaj [3,4]. A review of the literature on the fluid dynamics of solid particles in 
multi-phase systems showed different and incomplete approaches to the problem of the motion of various 
size particles near the wall of a turbulent fluid. In the work of Soo [10] a specific expression for the ratio 
of particle diffusivity and eddy diffusivity was found as a function of particle diameter, density, Reynolds 
Number (based on fluid stream turbulent intensity) and turbulence microscales. This expression is based 
on the analysis of the “probability of encounter for finite-size particles.” The final steps for calculating the 
ratio of diffusivities are summarized below for completeness. 
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Assuming that the Reynolds Number based on relative velocity between the particle and its surrounding 
fluid is small enough so that particle drag is given by Stokes’ Law, the particle diffusivity, D p , can be defined 
by the equation: 


,2 i 2 

Dp _ A* < u > 



with 


A* = A L /< 


> 1/2 , 


( 1 ) 


( 2 ) 


where < t/ > 1 / 2 is the intensity of fluid stream turbulence, and A^, and A^ are the Lagrangian and 
Eulerian microscales of stream turbulence, respectively. Introducing the impulse response parameter, J, 


I = 


2 

F A* 


( 3 ) 


The parameter I is defined as the ratio of particle impulse response time to the time a fluid particle remains 
in a velocity correlated region. For spherical particles, this parameter can be written as, 


I = (y/ir / 18) < Nju > (f>p / Pi ) (2r p / \ L ) 


( 4 ) 


where 


<N Re > = (2r p /u) < / (5) 

is the particle Reynolds Number. Upon substitution of Eqns. 2 and 3 into Eqn. 1, the ratio of particle 
diffusivity to eddy diffusivity becomes, 



Within the context of the BLIMPJ algorithm using the algebraic eddy viscosity models, the following 
simplifications and assumptions were made. 

Soo [10] considered A^ to be the radius of a pipe containing the two-phase flow. In our case, A# was 
assumed to be equal to the local boundary layer height, 6 . Since the turbulence intensity is not calculated 
in BLIMPJ as it is done with a Jfe — e turbulence model, certain simplifications had to be made. The eddy 
contribution to the shear stress is 

—p u* v* — r — p — CO 

dy 

From Schlichting’s book, the shear stress is related to the turbulence intensity by a correlation coefficient, 
given by 


From a limited set of data, ip was approximated to be —0.45. Then, from Eqns. 7 and 8, 


p) ^ = [e r / 0.45(/* + pe)] 1 / 2 


( 9 ) 


where e is the kinematic eddy viscosity. Since the problem usually is to obtain the ratio in Eqn. 6 for a 
given particle loading, it is essential to obtain an equivalent size spherical particle at a given point in the 
boundary layer. Particle loading is given by 


^p __ Pp n p 
W f ~ Pf 

Particle mass flow in an elemental ring located at radius, r, of width, A y, 
m p = (2 *r) (Ay) (p^) (MF P ) , ^ 

This is also equal to mass in the equivalent width of the particle phase ring, a, where 
mp = ( 2*r ) (a) (pp) , ^ 


( 10 ) 


( 11 ) 


( 12 ) 


From Eqns. 11 and 12, 
a - n p V p Ay 


(13) 


This formulation must be considered approximate since, clearly, a is not the radius of an equivalent sphere. 
A more detailed analysis is necessary for accurately computing the diflusivity ratio in such a situation. 

The significance of various terms in Eqn. 6 is dealt with in Ref. [10]. For a fixed X^/X E and small 
particle impulse response, /, a solid particle “follows” the fluid motion perfectly, and its diflusivity is equal 
to the fluid diflusivity. However, for a fixed A^/A^, a solid particle with large I does not respond to the 
fluid motion; it tends to remain stationary and not diffuse. 

In order to correlate the diflusivity ratio with the eddy viscosity, one can assume constant turbulent 
Schmidt Number. Consequently, 



This was coded in BLIMP J. 
CONCEPT CHECKOUT 


314 


In order to check out the turbulence model coded for the presence of particles in the boundary layer, 
the OTV nozzle configuration (Fig. 1) was chosen. As an example, metallic Aluminum was chosen to be 


the particle material present in the boundary layer. The particle loading was chosen to be W p /Wf = 0.5 in 
the whole boundary layer and BLIMPJ was run by turning on the new particle option. The boundary layer 
thicknesses for these runs including a clean flow run are plotted in Fig. 2. The A1 content in the boundary 
layer alone reduces the boundary layer height because of a reduction of total enthalpy of the system. The 
damping effects of particles in the fluid further reduce the boundary layer height because of a reduction 
in turbulence intensity. Profiles of EPSA (= p 2 e/p c /i c ) at a typical station 41 with and without particle 
loading are plotted in Fig. 3, where EPSA is seen to reduce in most of the boundary layer except in the 
middle. This is attributed to the species composition effect. However, the specific effect of damping only 
is to reduce the eddy viscosity. This has been observed by the authors in a separate calculation for this 
station. Profiles of the particle diffusivity ratios at three stations (13, 41, 48) are plotted in Fig. 4, where 
the effects are maximum at the throat. Comparisons were made for heat transfer rate in Fig. 5, where, 
again, the maximum change occurs at the nozzle throat. More credibility should be placed on this option, 
since it considers all the effects such as mixture, damping and “history” effects. In contrast, the previous 
option coded and described in Ref. [4] assumes constant property flow in the boundary layer, does not 
consider damping effects and does only “point” calculations. This formulation, however, does not include 
the computations for particle trajectory and, consequently, particle loading from station to station. 


THRUST DECREMENT CALCULATION 


BACKGROUND 

The thrust loss calculation originally implemented in BLIMPJ [1] has been modified by Praharaj in 
Refs. [4,5]. The modified procedure is particularly applicable to the projected OTV engine nozzles with high 
area ratio. 

Some of the assumptions made in the iteration procedure include: 


1. The average in viscid pressure obtained from TDK across the boundary layer width is a reasonable 
value to be used in BLIMPJ to define the boundary layer in the successive iterations. 

2. The longitudinal curvature effect is only approximated through the pressure averaging procedure. 

In the previous work [4], all these assumptions along with other standard boundary layer assumptions 
were made. It also provided procedures for calculating the performance of a rocket nozzle experiencing thick 
boundary layers for the two different cases given below: 

Case 1 - The potential nozzle contour is given and the objective is to define the hardware wall contour. 

Case 2 - The hardware wall contour is given and the objective is to define the potential contour and 
calculate the nozzle performance. 

The method had been checked out for the OTV nozzle using the procedure for Case 1. It used RAMP [11] 
as the inviscid code to define the pressure profiles and some of the pressure averaging was performed by 
hand. The current procedure uses TDK instead of RAMP to define the inviscid flowfield and a code (N WP 
code - Ref. [12]) was developed consisting of FORTRAN and control language to couple TDK and BLIMPJ 
to perform the iteration procedure in an automatic mode. The details of this code are given in Ref. [5] and 
will not be repeated here. 

In the above calculations the boundary layer thickness was measured from the wall in a direction normal 
to the nozzle axis. Since it is more accurate to measure the boundary layer normal to the nozzle wall, NWP 
was modified and the first iteration was completed. The zeroth iteration in BLIMPJ utilized the inviscid 
wall pressures from TDK and calculated boundary layer quantities including the thrust loss. This boundary 
layer width was sufficient to show a variation of pressure across the boundary layer region of the inviscid 
flowfield. Thus, an iteration was necessary in BLIMPJ to use an average pressure across the boundary layer 
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region obtained from the streamline information output by TDK. Note that the nozzle had to be extended 
at the exit plane using the exit cone angle in order to provide the streamline information in the vicinity 
of the exit plane. This is the first iteration. As an example, inviscid pressure profiles between the nozzle 
wall and boundary layer edge are given in Fig. 6. Iterations continued until a convergence of the thrust loss 
was achieved with a specified tolerance between successive iterations. Figure 7, on the other hand, gives 
the overall thrust loss of the nozzle as a function of the BLIMPJ iterations and compares the results of the 
new method described here against the more approximate method reported earlier. Method II, which is 
more exact than Method I, shows somewhat more thrust loss, but quicker convergence. The thrust vs. area 
ratio distribution from ODE, ODK, TDK and actual calculations using this procedure are given in Fig. 8. 
The differences between ODE and ODK are due to the kinetic losses. The differences between ODK and 
TDK are attributed to the two-dimensional losses. The thrust decrements due to the boundary layer losses 
were considered in generating the lower curve. Note that the thrust decrement iteration procedure described 
earlier has been taken into account for the lower curve. 

THRUST OPTIMIZATION 

The concept of using a high area ratio nozzle in the projected OTV engine was chosen because the 
specific impulse loss is minimized in such nozzles. As the area ratio grows in magnitude, two-dimensional 
losses become smaller and smaller. However, as the area ratio increases, so does the friction loss. This 
phenomenon has been pointed out by numerous investigators. As a result, the area at which the actual 
thrust reaches a maximum may be different from the exit area ratio. 

The subject of thrust loss optimization must be considered during the design stage. The desirability 
of improving on the “perfect” nozzle has been the subject of various optimization procedures. One of 
those is due to Rao [13] where, given the throat flow characteristics, nozzle length and ambient pressure, 
the procedure provides a nozzle contour yielding maximum thrust. This optimization procedure does not 
include the effects of wall friction. It was shown in Fig. 8 that as the area ratio increases, so does the 
boundary layer loss. Thus, the area ratio must be optimized to yield maximum performance. A thrust loss 
optimization procedure, which must be considered preliminary at this stage, but which includes wall friction 
without addressing issues such as structural design and other nozzle losses, is suggested below: 

1. Calculate a series of potential nozzles using such procedures as implemented in the Rao algorithm 
(or others) for various lengths. Each nozzle has a different area ratio and an optimized contour. 

2. Calculate the corresponding inviscid thrust for each nozzle. 

3. Once the nozzle contours are available, thick boundary layer calculations based on the procedure 
described earlier should be followed to define the thrust loss, and final displacement thickness distri- 
bution along the nozzle wall and the corresponding hardware nozzle wall contour. This defines the 
thrust loss variation with area ratio for nozzles with optimum contours. 

4. Optimize the resulting thrust ( = T inviscid — friction drag) with respect to area ratio subject to 
constraints such as length and structural weight. 

A parametric study was conducted for a typical high expansion area ratio nozzle such as the OTV 
nozzle, extended to an appropriate length beyond that necessary for the normal nozzle with exit area ratio 
of 1852 to provide the streamline information. The only parameter chosen for the study was area ratio. 
The geometry of the nozzle for each area ratio was assumed to be a part of the same extended high area 
ratio nozzle. The above optimization procedure was followed and the real thrust distribution along with the 
inviscid and boundary layer loss values are given in Fig. 9. 

In order to ascertain the ineffectiveness of larger area ratios (beyond a certain range) in producing 
thrust, Fig. 10 was prepared showing change of thrust from station to station as a function of area ratio. It 
is clear that beyond A/A* = 1852.22 where the change of thrust from its previous station was about 15 lbs., 
no more appreciable thrust gain was realized by further extending the nozzle. Although the original intent 
was to locate the zero for this curve, it was not achieved in this analysis for the following reasons: 
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1. It was not possible to run beyond A/A* = 1852 at REMTECH’s VAX because of space limitations. 

2. Since the in viscid thrust increases by extending the nozzle in a conical fashion, and the thrust loss 
also increases because of higher displacement effects with higher running lengths, the total thrust 
may be asymptotic in nature. This must be considered a conjecture at this point. 

3. On the other hand, if the nozzle were to be extended in a cylindrical fashion, the in viscid thrust 
should level off and the viscous losses should increase, thus resulting in negative change of thrust 
from station to station. This has not been exercised at this point. 

Thus, it is recommended that a separate analysis be made to define precisely the area ratio beyond 
which the thrust gain is negative. 


RELAMINARIZATION 


BACKGROUND 

The prediction of relaminarization phenomena remains as one of the strongest tests of validity of the 
turbulence models existing in the literature. Relaminarization is basically a reversion from turbulence to 
laminar boundary layer, principally caused by severe flow acceleration effects. Attempts have been made 
by various investigators to use the k — e turbulence model in an accelerating flow. These attempts achieved 
reasonable correlation with data. However, the thrust of the current work has been to develop an engineering 
model to update the existing turbulence models in BLIMPJ. Attempts were made in Ref. [3] to utilize the 
experimental data of Nash- Webber, which is one of the best-documented experimental investigations of 
compressible boundary layer. However, since this criterion worked only approximately for the Back and 
Cuffel nozzle data as shown in Ref. [3], the original formulation of Nash- Webber was closely examined. The 
following modifications were made (Ref. [5]) for including the effects of wall cooling and roughness in the 
transition criterion curve-fit. 

1. The Nash- Webber correlation was valid for adiabatic wall condition. 

2. The transition line was moved down for cooled walls, i.e., it is easier for turbulent flow to laminarize 
on cooled walls. 

3. The transition line was moved up for rough walls, i.e., the phenomenon is the opposite of cooled wall 
flow. 

One other important factor that was not considered in the above study was the effects of free-stream 
turbulence in the boundary layer which would oppose relaminarization. This has been documented in Ref. [5] 
and is summarized here. 

IMPACT OF FREE-STREAM TURBULENCE 

Free-stream turbulence in liquid rocket nozzles is caused by the violent mixing at the injector and by the 
explosive burning of the fuel-oxidizer mixture in the combustion chamber. The role of free-stream turbulence 
in gas-turbine systems has long been recognized and has been measured in certain situations. These studies 
indicate that the primary effect of an increase in free-stream turbulence is the upstream movement of the 
onset of transition. The effects of free-stream turbulence on turbulent boundary layer profiles have been 
found by others to give slightly fuller profiles and higher turbulence levels, resulting in higher momentum 
thicknesses, smaller form parameters, and increased heat transfer as well as skin friction coefficients. In the 
rocket nozzle situations, where relaminarization is a possibility because of high favorable pressure gradients, 
free-stream turbulence will delay relaminarization. 

In the current BLIMPJ framework, a lot of simplifications had to be made to integrate the effects 
of free-stream turbulence into the algorithm. The relationship between turbulence intensity and Reynolds 


shear stress has been dealt with in various references including Hodge and Adams [14]. The proportionality 
between these two quantities depends on whether one is considering the inner or the outer layer and also on 
the nature of the boundary layer, incompressible or compressible. Roughly speaking, for isotropic turbulence, 


u'v' = y> u ' 2 ( 16 j 

where the correlation coefficient, V, is approximated to be -0.45. The free-stream turbulence level, T u , is 
given by, 



Thus, from Eqns. 16 and 17, 

) =xbU?T 2 /iq\ 

V / Boundary Layer Edge ™ e u (1®) 

This was converted to an equivalent EPSA (= p 2 e/PeM«) value in subroutine TRMBL with the proper non- 
dimensionalization at the boundary layer edge for each station, and was ramped down to a value of zero at 
the beginning of the outer layer of the two-layer eddy viscosity model. The OTV nozzle was used as a test 
case for a constant value of U e T u = 200. Since U e increases continuously in the nozzle as the flow expands, 
the value of T u drops off from 107.0 percent at the inlet to 1.24 percent at the nozzle exit, as shown in 
Fig- The variation of EPSA is given for station 6 in Fig. 12. A plot of heat-transfer rate was made in 
Fig. 13a for the above variation of free-stream turbulence giving slightly higher heat transfer value at the 
nozzle throat than the normal quiescent situation, as seen from Fig. 13b. The drop in heat transfer ahead 
of the throat in the nozzle convergent section is not clear. 

The approach described above, even though innovative, must be checked out against other codes and 
any available data before it can be used as a reliable engineering tool. 


RECOMMENDATIONS 


The work presented here represents the last part of the work performed under NASA Contract NAS8- 
36551 from Marshall Space Flight Center, Alabama. The following are a set of recommendations in support 
of the present work. 

1. The various modules coded in BLIMPJ should be checked against other available data. 

2. Since full-fledged CFD codes are becoming more and more a reality, the ideas in the above modules 
should be included in the CFD codes and checked against the JANNAF codes and experimental 
data. 

3. New tests should be devised and conducted to provide a better database in the areas of particle 
effects and relaminarization studies. 

4. More test data are necessary for validating the wall roughness effects and thick boundary layer 
studies. 
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Figure 1: OTV Nozzle Configuration and Turbulent Boundary Thickness (Selected Stations 
Shown are for the Condensed Phase Study). 
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Figure 3: Profiles (typical) of EPSA with and without Particle Loading at Sta. 41. 
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Figure 6: Inviscid Pressure Profiles Between Nozzle Wall and Boundary Layer Edge in First Iteration. 
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Figure 12: EPSA Profile for a Station of the Convergence Section. 
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